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ANALYSIS OF A FINITE DIFFERENCE GRID 


Goetz H. Klopfer 

Nielsen Engineering & Research, Inc. 

SUMMARY 

Some means of assessing the suitability, of a mesh network 
for a finite difference calculation are investigated in this 
study. This has been done by a study of the nonlinear trunca- 
tion errors of the scheme. It turns out that the mesh can 
not be properly assessed a priori. The effect of the mesh on 
the numerical solution depends on several factors including 
the mesh itself, the numerical algorithm, and the solution. 
Several recommendations are made with regard to generating the 
mesh and to assessing its suitability for a particular numerical 
calculation. 


INTRODUCTION 

One of the most important problems that arises in finite 
difference solutions to physical problems is the- quality of the 
grid used in the calculation. For simple configurations, e.g. , 
two dimensional problems, a conformal transformation which gives 
an orthogonal grid can be used. This still leaves the questions 
of the adequacy of the grid clustering unanswered. In complex 
problems the grid may not be orthogonal, and a question of the 
effect of grid skewness on the solution accuracy arise. In 
addition to these problems of grid clustering and skewness, it 
is desirable to determine whether the grid is sufficiently fine 
to capture the detailed physics of the flow, especially if the 
location of phenomena like shock waves change during the solu- 
tion. Other critiera (e.g., smoothness) may also be important. 

At present the only way of determining the suitability of 
a grid is by visual examination, and while this may be 
satisfactory for simple geometries, it is an almost impossible 


task for complex geometries. It would be extremely useful, 
therefore, to derive some criteria to measure the suitability 
of a grid that could be determined computationally. The 
quality of a complex three dimensional grid could thus be 
determined before a lengthy flow solution is attempted. 

In a general two dimensional grid the coordinate trans- 
formation T “ t, C = C(t,x,y), n “ n(t,x,y) permits clustering 
and spreading of the points in two-dimensions. The components 
of the Jacobian transformation matrix (hereafter called 
"metrics") Xj.,y^ and x^,yj, can be examined to determine the 
quality of a grid. The main problem is to decide the attributes 
of a good grid; it is anticipated that quantities such as 
curvature variation, Jacobian variation, skewness and volume 
of grid cells, in addition to other geometric quantities will 
be important. 

First the attributes of a good grid are determined, and 
second, means of quantifying the grid quality are developed. 

In addition to this purely geometric problem, it is likely 
that the relative quality of a grid will depend on the algorithm 
in the code since the major aim is to reduce the truncation 
error. Since different algorithms have different truncation 
errors it is probable that grid quality is dependent on the 
algorithm and the solution. A particular algorithm is chosen 
for more detailed study. 

Finally, a feasibility analysis on the extension of the 
above criteria to three dimensions is undertaken. 

COMMENTS ON MESH CRITERIA 

Before some criteria can be established which will serve 
to quantify a "good" mesh, it is necessary to discuss the various 
purposes a mesh needs to serve. This will be done in this 
section. The discussion will be limited to problems applicable 
in computational fluid dynamics. 
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For an ideal incompressible flow, finding the solution is 
tantamount to finding the coordinate lines. This is of course 
the motivation of conformal mappings. The solution and the 
mesh are given in terms of level curves of potential and stream 
functions. Extensions of this idea to compressible flows has 
not been as successful. For example, the "contour dynamics" 
scheme of Harlow (Ref. 1) recast the problem as finding the 
locations of level curves of density, pressure, and velocities. 
While the problem of establishing the mesh is circumvented, 
other more severe and less tractable difficulties arose. 

Another method of avoiding the mesh generation problem is 
by the use of Lagrangian variables (Ref. 2) . In these 
variables the mesh node points move with the fluid velocity 
and thus the mesh is automatically generated. The unsolved 
problem with this approach is the excessive mesh distortion 
at stagnation points or in the viscous layer. This distortion 
causes accuracy and stability problems. Attempts have been made 
to solve this problem by remeshing before loss of accuracy 
occurs. These methods are called Euler-Lagrange methods by 
the Los Alamos school (Ref. 3) . 

Other uses of meshes for numerical algorithms are to reduce 
the trunaction errors of finite difference schemes. In other 
words, the mesh of fixed number of node points is adapted so as 
to improve the accuracy of the numerical solution. This has 
been the motivation of the work of Brackbill (Ref. 4) , Klopfer 
(Ref. 5) , and others. Hindman (Ref. 6) established certain 
properties that the transformation and the metrics must sat- 
isfy for certain mesh dependent truncation errors to vanish. He 
showed how these differ if the governing equations are written 
in various different forms, e.g., strongly and weakly conserva- 
tive forms to name only two. This indicates that the form of 
the governing equation can also have an effect on establishing 
mesh criteria. 


3 



















If conservation is important, as is usually the case for gas 
dynamics, where shock waves and tangential discontinuity surfaces 
are to be captured by the numerical scheme, then the mesh must be 
such that numerical conservation is maintained to a certain level 
of approximation. This requirement precludes rapid changes in 
mesh cell size. The changes must be smooth and gradual if con- 
servation is required. 

Another use of adaptive meshes is to keep track of fluid 
interfaces such flaime fronts Dwyer (Ref. 7) and water waves. Here 
the accuracy of the solution is strongly dependent on the proper 
resolution of the flame front temperature gradient or air-water 
interface. The truncation error of the flow variables is not the 
problem here. The difficulty is obtaining the source terms 
accurately for the chemical kinetic equations or the boundary 
surface for the wave problem. The mesh effect on the accuracy of 
these flow fields will be different than those based purely on 
truncation errors. 

Another purpose of a finite difference mesh is to control 
the stability, well-posedness, or convergence properties of num- 
erical schemes. For example, Hagin (Ref. 8) used this method to 
keep an integral equation method well posed. Preconditioning 
methods or multi-grid methods also make indirect use of this 
mesh property to increase the convergence rate of numerical 
schemes Lomax (Ref. 9) . 

From this brief survey it is obvious that there are too many 
different requirements that a mesh has to satisfy. To establish 
criteria based on all these functions of the mesh is a far too 
ambitious undertaking. For this reason the remainder of this 
study will focus attention on meshes adapted to reduce the trun- 
cation errors only. The governing equations considered will be 
in strongly conservative form with no source terms due to physi- 
cal processes. In the transformed or computational space the 
equations considered will be both strongly and weakly conserva- 
tive. The source term for the latter form will be due to the 
transformation metrics only. 
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In the next sections the treinsf orrnation , governing equations 
and truncation errors of a simple numerical scheme are derived. 


CURVILINEAR COORDINATE SYSTEMS AND 
GOVERNING DIFFERENTIAL EQUATIONS 


Metric Tensor 

In any two dimensional surface with local coordinates 
5^ = CC*nl we have a local element of length 


ds^ = g^.d^^dC^ i,j = 1,2 

= + 2gj^2'^Cdn + g22*^n^ 


( 1 ) 


If imbedded in -a plane "physical" space with coordinates 
x^ = (x,y) then the components of the metric tensor g^^^j are 
given by 


^11 = 


^12 




( 2 ) 


2 2 

g_. = X + y 
^22 n ■'n 


so that ds 
is 


2 2 

dx + dy . The determinant of the metric tensor 


g = Ig. .1 = j2 = g^^g ^2 - 


(3) 


where J is the Jacobian of the transformation between the two 

coordinate systems. If the map; ing between the two systems is to 

-1 

be one-to-one, then J 0 and J / 0 must be satisfied. 
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With this transformation a mesh can be laid out with the 
node points at (.convenient) integer values of (5/n)- This mesh 
is simply a cartesian mesh with square cells and will be taken 
to be the computational mesh and (5»n) are the computational vari- 
ables. The metric tensor components can be given a physican 
interpretation. The g^j can be rewritten as 

= J^(Vn • 7ril 

gi2 = • 70 (4) 

922 “ * 

Thus g^^j^ gives a measure of the length (squared) of the 
side of the cell parallel to the 5 coordinate, 922 the 
length of the side parallel to the n coordinate, and g ^^2 is a 
measure of the orthogonality between the cell sides. The metric 
tensor components however do not completely describe the trans- 
formation between the two spaces. The orientation of the (C,n) 
coordinate system relative to the (x,y) system is not fixed by 
g^j and must be specified separately. One such specification 
can be the angle between C-axis and the x-axis as given by 

cos<> = ^ (5) 



Other representations are possible but this one is simple 
and with a surface conforming curvilinear coordinate system it 
gives the angle between the tangent plane of the surface and the 
x-axis. It will be called the "solid-body" rotation angle. 

This solid body rotation is only defined for flat 2D surfaces. 

On a curved surface the integral of the rotation about a closed 
contour is equal to the integral of the curvature of the surface 
enclosed. 

I 
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For later reference is useful to have the transformation 
metrics in terms of g. . and (P. 


and 





COS(tl 

sin>p 

(gi2 cos 4 > + /g sin<Ji) 
(/g cos!j( - gj^2 sin({i) 


( 6 ) 


The latter two become for an orthogonal system * 


^ = "’^22 
Yjy = /?22 cosifr 


Note the strong dependence of the metrics on the solid body 
rotation. 

The metric tensor can be rewritten or split in various 
different ways. For example, 




a coshB sinhS 

sinhS — coshB 
a 


where a is the cell aspect ratio, B is a measure of orthogo- 
nality and J the cell volume (area) 


and 


ct = >^gii/g22 
sinh = 


7 
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Or can be split into dilation and shear components 

tTi 0| -fsinG COS0 "] 

^ij 2[0 ij ^[cos8 -sinej 

where the first matrix represents a viniform dilation and the 
second the shear. T and R are the strength of the dilation and 
shear respectively. T is, of course, the trace of g.., that is 


^=^ 11 ^ ^22 


also 



r = 

and 

sine = A/7 a^ + B^ 

where 

A — 

and 

B = gi2 . 


ORIGINAL 
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With these definitions of the metric tensor, certain mesh 
criteria can be established. 

MESH CRITERIA 

In recent years curvilinear coordinate systems have become 
important for improving the efficiency of numerical solutions of 
flow fields around arbitrarily shaped bodies. In generating 
these meshes many workers have incorporated various properties 
they considered desirable into the mesh. It is not clear why 
these properties are required for a mesh to be "good" for a 
certain numerical scheme. 

Brackbill (Ref. 4) generates two dimensional meshes 
such that several of the mesh pro^-i. ties are simultaneously 
optimized. For example, he considers cell volume, J; the trace 
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of g. j which he calls cell smoothness, T; and orthogonality, g ^2 
to be important. In other words he generates a mesh such that 
the conflicts among the three requirements are minimized. 

Yanenko (Ref. 10) also considers orthogonality and a specified 
variation of cell volume to be important. However, he does not 
use the cell smoothness, instead he considers a mesh "non- 
Lagrangian-ness” which is the difference between the fluid veloc- 
ity and the mesh node velocities. Steger, et al, (Ref. 11) have 
developed a very rapid mesh generator by satisfyi.'\g only two 
properties, for example, orthogonality and a simple variation 
of cell volume. Ablcw (Ref. 12) and VThite (Ref. 13) use the 
solution in a direct manner to generate the mesh. In other words, 
the solution and mesh are solved simultaneously such that the 
measure of "difficulty” of the method is spread out uniformly 
over the computational domain. The "difficulty" can take on 
several meanings, for example, equal arc length along the solu- 
tion curve, truncation error, single step error, or a number of 
other things. If the differential equation is the Poisson 
equation and one takes the truncation error to be the measure of 
"difficulty", then Ablow (Ref. 12) consider 5]^^' ^12' ^22' 
along with the higher order derivatives of .he solution to be 
the important mesh criteria. Thompson (Ref. 14) have considered 
arc lengths between coordinate lines (i.e., g^ , and ?22^ 
orthogonality as mesh criteria. These six works are probably 
the state-of-the-art in mesh generation with certain desirable 
properties explicitly built into them. 

It should be apparent that if the mesh is solution adapted 
and the solution is a vector, then there should be a mesh for 
each component of the solution vector. This of course involves 
the additional expense of generating and storing multiple meshes 
in addition to the extra burden of interpolating eimong the 
meshes. These types of difficulties were experienced by Harlow 
(Ref. 1) with his "contour dynamics" method. .While multiple 
meshes may be unnecessarily cumbersome and expensive they should 
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not be rejected out of hand for they may be advantageous for 
certain kinds of problems. But, in general, a single solution 
adapted mesh is generated based on some scalar composite of the 
solution vector, e.g., pressure. 

It is not evident that any of these properties are in any 
way desirable for a "good" mesh. The very notion of a good mesh 
is vague. The function of any mesh network, dynamic or station- 
ary, is to discretize a system of governing partial (usually, 
non-linear) differential equations so that they can oe solved 
by a numerical algorithm. In this report only finite difference 
schemes will be considered. For finite difference schemes the 
m.esh should minimize the truncation error for a given n'lmber of 
discretization points. This is because one can not easily add 
or delete mesh points during a calculation. Since the trunca- 
tion error depends on the transformation metrics, the solution, 
and the numerical scheme, it may be expected that the desirable 
mesh properties will be scheme and solution dependent. 

It is also possible that the mesh criteria may depend on the 
order of the differential equations. For this reason twc types 
of differential equations will be considered. The first type 
will be a system of nonlinear first order partial differential 
equations, such as the Euler equations of gas dynamics. The 
second type will be a scalar second order partial differential 
equation such as the steady full potential equations for tran- 
sonic flows. 


TRUNCATION ERROR 

First Order Differential System - 
First Order Scheme 

In the following sections the truncation error of two 
algorithms on an arbitrary curvilinear mesh will be derived. 

From the truncation error the desirable form of the transforma- 
tion metrics will be extracted. The governing equations will be 
conservation law form 
10 
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where f and g are nonlinear functions of w. 

A generalized coordinate transformation is introduced 

T =» t 

: C = Ct,x,yl (8) 

n = (t,x,y) 

The transformed governing equations are written in two forms, 
the strong conservation law (SCI,} and the wea)c conservation law 
(WCL) forms: 

Strong conservation lax: (S CL ) 




(wj) = - xw) - x^(g - ^)j 

(f - xw) + (g - a 0 


(9) 


Weak conservation lav (WCL) 


riw 0* -"1 f/* ♦ — — — ) 

TF * ^ 

+ ^ - yx^)w - y^f + x^?| 


g = 0 


( 10 ) 
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To both of these forms the mesh velocities (x,y) must be given 


- X = 0 

“ y = 0 


( 11 ) 


For convenience in determining the truncation error the mesh 
accelerations (x,y) are also given 


(x) ^ - X = 0 

(y) ^ - y = 0 


( 12 ) 


Here x, y, x, y are given functions of w, w , w , w,y. 

Both of these forms are written in terms of new dependent 
variables and fluxes as 


\ 


I 

I 


”t * *5 * \ 


+ h = 0 


( 13 ) 


For the SCL form 


w = 


g = 


WJ 

X 

y 

« 

X 

y 


-y,(f 


f = 


y^ (f - xw) - x^ (g - yw) 
0 
0 
0 
0 


xw) + (g - yw) 
0 
0 
0 
0 


h = 


(I'D 


12 
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and for the WCL form 






xw) + x^(g 


(last three terms of equation (10) 


Solving system (13) by a first order mmerical scheme 


,n _n 


„n+l . „n f -__fi ^ gj^-1 ~ ..ll + = 0 

Ir 


one obtains the following modified equation 


+ h + ^ ^ £25 + ^2"^ 92 n = ^ 


= w . This notation will be used throughout this 
TT 


where w^^ 
report. 
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By eliminating in terms of spatial derivatives one 
obtains 


''t * h ^ ^ ^ (-S,) ■ aah 


2 2 2 2 
where 0 (A ) represents 0(Ac.,A ,A,,Aj.A ). 

t n T ^ n 

Here s = + h and (-s^) is given by 

C-s^) = [h^ + + ^ww/cc ^ww ''nc W ^ww/5n 

i n n 4 

+ c w [•s+lf +h +f w_ + f_ w_c 
'ww^ nnj L w w^w C 2w^ 25 

+ f w _ + g w + g_ w_ + g w_ I • s^ 

w^w^ n5 ^w^w n ^2w^ Cn 5 

+|g +h +f w-+f w_^ + w . + a w 

r=w w w w S w w_ . 2w ri5 -w w n 
I- n n n5 n n 

+ a w_ + o_ w- • s + Ih + f I • Spp 
'W^W^ 5n -2w^ 2nJ n [ wj 55 

+ fh +f +gl*S[.+rh +ql*s 

The meaning of these terms should be made clear. Terms with 
subscripts 5 or n» e.g., s^^ , are 


s. = rV-(s) = + g., + h) . 

^5^ ^5% ^ 


Terms with subscriots w, w_, w , etc. are the Jacobian coeffi- 

5 n 

cient matrices or flux Jacobians (not to be confused with the 
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Jacobian of the coordinate transformation) . For example 


ww. 


8^f 

>w3w, 


This term is a third order tensor containing 8x8x8 terms if w is 
a 4-component vector as for example in the two-dimentional 
equations of gas dynamics. 

As can be seen this modified equation is a rather complicated 
function of the spatial derivatives of the metrics and the flow 
variables. It is far too complicated for one to readily extract 
any useful information. It is presented, however, to emphasize 
the additional effects that a dynamic mesh has on the truncation 
error. Thus one should not be surprised by difficulties experi- 
enced with dynamic meshes, as used for example by the Space 
Shuttle codes Kutler (Ref. 15) or the parabolized Navier-Stokes 
codes, Schiff (Ref. 19). 

If we specialize the problem to an unsteady flow on a fixed 
mesh then we obtain 




|h + f w- + g w I 
2 ( w ww 5 ww r,j 


S + f„s^ + g, 


w®nj 


O(A^) 


(19) 


Note that for this case all of the coefficient matrix Jacobians 
derived from the spatial derivatives of w vanish. 

Comparing the two modified equations immediately reveals 
that the truncation error contains derivatives of one higher order 
for the dyncunic mesh. Thus determining proper mesh criteria are 
more restrictive for a dynamic mesh than for a fixed mesh. If we 
further simplify the problem to steady flow then we obtain 


15 


f^ + g^ + h + ^£2? - ^2^1 “ <20) 

There is no further reduction in the order of the derivative 
in the truncation error. The equations are recast as 

[**^^Je* 

White CRef. 13 } speaJcs of an ideal mesh for a one-dimensional 
problem as one where the numerical and exact solutions are iden- 
tical. That is the truncation error vanishes. The ideal mesh 
for the above first order numerical scheme is one that eliminates 
both of the truncation errors. That is, the mesh is such that f 
and g are linear in 5 and n, respectively; 

f = + K2 

g = K3H + K4 

where = -K^ = constant and K2 and are arbitrary functions 
of C and n, respectively. 

For the SCL Form 


f = yr^f - x^g 
g = -y^f + x^g 

Thus a system of differential equation for x and y in terms of 
"known" solutions T and g is obtained. 

* ’‘2"' 
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This, of course, has no solutions for nontrivial 1 and g. Thus 
the ideal mesh does not exist. The effect of the metrics on the 
trtincation error can be seen by examining the trioncation in ex- 


panded form. For the SCL form the errors are ^nd ^ 


2^nn' 




^^55 - \Hk 


Similarly 


= -y_ J + Xr g + 2{-yp 1 + x- g )-yr^ + Xj-o^ 


The truncation errors for the SCL form depend linearly on 
the metrics and their derivatives. The errors do not depend in 
any direct way on the Jacobian, orthogonality or cell smoothness 
parameters except through the metrics. They do, however, depend 
more directly on the solid body rotation and cell side lengths, 
g^^^ and Q 22 (individually, not combined into a smoothness 
parameter) . 

For the V7CL form the errors are 


• ■ (Sa - (W A • &K - (Si.. 
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and 


=nn 


(^‘l’ • (?b • ■!(% • (%1 

• (^>, • (?K 


(25) 


The errors here depend linearly on the metrics and inversely 
on the transformation Jacobian and their derivatives. It can be 
expected that the WCL forms will require meshes with smoother 
Jacobians than the SCL forms of the governing equations. This 
conclusion, of course, assumes that the metrics in the source 
term are computed exactly. 

For both forms the mesh cell aspect ratio does not enter into 
the individual tr\incation errors. However, if the error along one 

coordinate line is much larger than along the other then the ratios 
Yc X. 

~ and ^ need not be near unity. The ratios of the metrics for 

an orthogonal mesh are determined by the cell aspect ratio and the 
solid body rotation angle. 



tan<!» 


^22 


Fixed mesh-unsteadv flow 

In this section the expanded form of equation (19) will be 
examined. To keep the algebra within reasonable bounds, the 
governing equations will be limited to the linear scalar convec- 
tion equation, i.e., 


u^ + c u + c u 
t XX y y 


0 
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where c and c are constant convection velocities in the x- and 
X y 

y-direction, respectively. For this exeunple the over barred 
quantities corresponding to equation (7) are 


w = u 


£ - c^u 

? = Cy» 


The transformed equation is given by equation C13) and for 
a fixed mesh the transformed variables are for the SCL form 

SCL 

w = Ju • 


( c y c X \ 


Ju 


and 


/"C„v, c.x.- 

- (-^ * 


g 

h = 0 


The flux Jacobian matrices become 


, _ , 
w J 




= V 


(26a) 


and 


fww = ° ' 


g = 0 . 


Therefore the modified equation is now 


w + 


+ ^f 

2 ^2, 


in. 


At U 


^ T^2 ^ — pw ' ^ 

r» ' ' 


I = 


nl 


C(i^) 


(26b) 
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where 
and 

* «2n 

In expanded form the modified equation is 
Ju^ + (JX)u^ + (Jv)Uj^ + * 2(J^)5U^ + 

+ ^|(JV)U2^ + 2(Jv)^u^ + (Jv)2^u| H- ^|x[(JX)U25 (JX)^u^ 

+ (JX) (v) [2u^j| = O(A^) . (27) 

where X and v are given by equation (26a) . 

There are several ways this equation can be simplified for 
particular cases. Suppose the mesh is to be such that it induces 
no truncation errors when the flow is uniform. This is called 
the uniform free streaun test by Hindman (Ref. 6). Applying the 
uniform free stream test u - 1 yields; 

- Vn’2c| ^ " V5^2n|= ^^8) 

This equation gives the first order triuication error due to 

the mesh above for a uniform solution u = 1. !Jote that there is 

no mesh error with a At dependence as all those terms cancel. 

2 

This equation also gives the requirement that u = 0 (A ) , that is 
the metrics must be smooth enough such that 

" Vn^2, ^ 
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and 

'"Vc * Vs’2n ' 

If these conditions are not met then the free stream value 
of u = 1 will not be recovered to the formal order of accuracy 
of the numerical scheme. This uniform flow test is one method 
of testing the adequacy of a discretization mesh and is similar 
to the approach taken by Hindman (Ref. 6). Hindman's approach 
differs in that he tried to satisfy u^ = 0 by special differ- 
encing of the metrics for the various form of the governing 
equations. In the present approach the metrics are differenced 

no differently than the flow variables and it was only required 
2 

that ) for the first order numerical scheme used. 

It is quite obvious that the uniform free stream test is 
artificial and somewhat meaningless for a nonuniform flow field. 
If the flow field is uniform there is no need for the numerical 
scheme. As shown by equation (27) the mesh dependency of the 
truncation error appears in a much more complicated manner and 
it is somewhat difficult to separate the mesh truncation errors 
from those caused by th^ -olution. This will be attempted in a 
later section for a second order scheme and a steady state 
solution. 

For the WCL form of the transformed equation, the trans- 
formed variables are 

WCL 


w = u 



21 



origina!- pass « 

OF POOR QUALITY 


and 




The Jacobian raetrices axe now 


W “ J J “ 


Z£i.!z!li = 


»• h/u 


(u ?f 0) 


and 


g = g -In. =0 
^ww ^ww ww 


The modified equation is now 


w^ + + h + 


~T^2Z T^2ti 


A^|h - s + f *Sp + g •s| = 
2 ( w w 5 ^w n) 


O(A^) 


where s = + h. 

This equation may be expanded to obtain 
u^ + (X)u^ + (u)u^ + ^|(X)u.^ + 2(X^)u^ + (X)^ u| 


+ (Xv). 


( 29 ) 


•\! 


= O(A^) 


( 30 ) 
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For the WCL form the uniform free stream test gives 


“t ^ 2 


l^x^n ~ °v^nl 
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. I “C Vr + C Xr| 

= OCA) 

2 I -J '2n 


(31) 


This result is similar to the SCL form equation (28) . The 
requirements for no mesh error are now 


and 
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-v- 0(A) 


-°xy5 


2n 


-V 0(A) . 


These conditions are probably more restrictive than for the 
SCL forms in that the Jacobian must also be smooth. These con- 
clusions are similar to those reached by Hindman (Ref. 6) . Here 
we conclude that the optimum mesh (optimum based, e.g., on the 
uniform free stream test) depends on the form of the transformed 
equation used. Hindman concludes that to recover the free stream 
values the metrics must be differenced differently depending on 
the form of the transformed equation solved. Of course, the mesh 
based on this test is not optimum in the sense of minimum trxinca- 
tion error of the numerical solution. For this the flow field 
solution and the lower order derivatives of the metrics must also 
be considered. 

It is of interest to see how the truncation errors and thus 
possible mesh criteria change as a higher order numerical scheme 
is used to obtain the numerical solution. The derivation of the 
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modified equations both in general and expanded form for a second 
order accurate leap frog scheme and system tl3) have been rele- 
gated to the Appendix. The major effect of the higher order 
scheme can be seen by comparing the error terms of equation (A-5) 
to the errors given by equations (22) and (23) . The major dif- 
ferences, besides the additional terms, are the appearance of 
derivatives of one higher order for both the solution and the 
metrics. This indicates that for higher order schemes, the 
higher order derivative of the metrics must be smoother than for 
a lower order scheme. 

Dynamic mesh — unsteady flow 

The expanded form of the modified equation for the dynamic 
mesh will not be given. The algebra is nearly intractable. A 
complete study of the dynamic mesh case requires the use of 
machine algebra and that is beyond the scope of this study. It 
is also doubtful whether anything useful in way of a solution 
independent mesh criteria will evolve from such a study. It is 
clear that it is not possible for the stationary mesh case covered 
above. 
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Second Order Differential Equation — 

Second Order Scheme 

In this section a scalar second order partial differential 
equation solved by a simple second order accurate numerical 
scheme will be examined. The steady state version of the full 
potential equation for transonic flow is the model equation, 
Holst (Ref. 16) . 

The full potential equation written in strong conservation 
law form is 

(0<(>^) + (P<j>„) = 0 

^ X ^ y 
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where the density, p, and the velocity components $ , ^ , are 

X y 

nondimensionalizod with respect to the stagnation density, 
and the critical sound speed, a*, respectively; x and y are the 
Cartesian coordinates and y Is the specific heat ratio. 


In the transformed space (given by C • C(x»y) and n • n<x,y)) 
the equation becomes 


whore 


(JoU)^ + (JpV)^ - 0 
U - (922$^ - 

V - 


(33) 


The definitions of 9^-,* 9 .,-m and J are given by equa- 

tions vl') and \3). The density is now 



It should bo noted that for this second order differential 
equation that the "traditional" mesh criteria reappear directly, 
namely, J, 9j^^» 9 t>, and This is unlike a first order 

differential system where only the metrics and J appear. 

Since for full potential equations the density is usually 
the desired -'olution from which the pressure and Mach numbers 
can bo dotoi.«iined we will examine the truncation error of 
equation (34) only. For simplicity simple central differencing 
will be used throughout the entire flow field whether the local 
.Mach numbers are sub- or supersonic. Let the "exact" density bo 
given by Og and tliO numerically derived density by . The 
relation between the two is 
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pg ■ pj^+ Ap + higher order terms 



The Ap can be derived from the modified equation form of 
equation (34). Replacing all the derivatives by central differ- 
ences and expanding in a Taylor series obtains the modified 
equation of (34) as 



reremi i, r „ ^ r „ i . T^nTf n. 
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.(S • ■ SJ.. 

and 

4J - - x^yj^) . 

Note that the effect of the trvmcation errors cf.the metrics 
and the solution (<>) are somewhat separable. In equation (35) 
the first group of terms on the right hand side is the truncation 
due to the solution. This part of the error still has the metrics 
as coefficients. The second group of terms is the mesh induced 
truncation error and also depends on the solution. Only in this 
sense are the two errors separable. Thus, for a given solution 
one could check the adequacy of a mesh by computing the two 
separate errors as above and require that the mesh induced errors 
be no larger than the solution induced errors over the entire 
solution domain. 

This is not a very satisfactory result in that a mesh in- 
duced truncation error only has meaning with respect to a parti- 
cular flow field solution. A numerical example for the full 
potential equation is presented in a later section. 

MESH- INDUCED NUMERICAL INSTABILITY 

The use of curvilinear mesh or a nonuniform Cartesian mash 
can introduce what will be called mesh-induced numerical insta- 
bility. Although a numerical scheme may be stable for a uniform 
Cartesian mesh, there is no guarantee it will remain stable on 
another type of mesh. This assertion can be quickly verified by 
o.xami."'ing the modified equation (A-5). In this equation there 
are four dissipative terms whose coefficients depend on the local 
metrics and their derivatives. These four terms are 
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These terms can act as either deimping or destabilizing terms 
depending on the signs of and If they are destabilizing, 

then they must be stabilized by either the inherent dissipation 
of the scheme (which the leap frog scheme does not have) or by 
the artificial dissipation usually appended .o t’^e scheme. In 
some cases for which the metrics vary rapidly enough, the num- 
erical or artificial dissipation may not be sufficient and thus 
a mesh induced instability can occur. An example of this occur- 
rence has been shown by Hindman (Ref. 6). 

This unpleasant problem does not occur for uniform Cartesian 
meshes and no solution is offered here. We simply wish to raise 
a warning flag that even though a mesh is generated according to 
some criteria based on truncation errors there is no guarantee 
that the numerical scheme will be stable. 

NUMERICAL RESULTS 

Some numerical results will be presented here. The results 
will be in terms of a post-mortem. For a given mesh and numeri- 
cal solution based on the full potential code (Holst, Ref. 16) 
as modified by Dougherty (Ref. 17) the truncation errors will be 
analyzed to see if the mesh indeed was adequate. 
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The grid is shovm in Figxire 1. Shown is a C-type mesh 
around a NACA 0012 airfoil. The far field boundaries (not shown) 
are approximately six chord lengths from the airfoil in all 
directions. It consists of 175 ^ 31 node points and is generated 
by use of the Poisson equation with some grid control (Sorenson, 
Ref. 18) . A solution in terms of density is shown in Figures 2 
and 3 for a free stream Mach number = 0.80 and angle of attack 
0=0®. These results are presented as carpet plots with the 
computational variables I, J (representing n, respectively) 
as the independent variables. The airfoil surface is given by 
J = 31 and the far field boundary by J = 1. The leading edge of 
the airfoil at I = 88 and the trailing edge at I = 31. The plots 
are distorted but this should cause no confusion. The viewpoint 
of Figure 2 is thus at the far field boundary and the leading 
edge. In other words, from the right hand side of Figure 1. The 
viewpoint of Figure 3 and all the subsequent figures is from the 
far field boundary and the trailing edge, or the lower left hand 
corner of Figure 1. The rise of the density at the leading edge 
stagnation point and at the shock wave is quite apparent. 

The next set of three figures presents some of the tradi- 
tional mesh criteria. These are the smoothness ^ 22 ^ 

Figure 4, the Jacobian (plotted inversely) in Figure 5, and the 
orthogonality (?]l 2' "skew") in Figure 6. The plots 

have all been normalized with the normalization constants shown 
on the plots. These figures show that the smoothness parameter 
is very uniform near the airfoil and becoming less so at the 
far field boundary. The Jacobian plot. Figure 5, indicates that 
most of the variation occurs near the leading edge and relatively 
little elsewhere, T)ie mesh is nearly orthogonal except at the 
trailing edge and far field boundary, as shown by Figure 6. 

Figures 7 and 8 present the variations of g^^^^ and 922* 
shown they are nearly uniform except at the trailing edge and 
far field boundary, wit.h g 22 being an order of magnitude more 
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• important (its normalization factor is 10 times larger) . The 
solid body rotation, cos<t, is shown in Figure 9 and is labelled 
as ROTAT. The cell aspect ratio, ^9n/922' ASPECT, is 

given in the following figure. Both of these indicate large 
variations near the trailing edge and the leading edge. 

The following set of four figures present the metrics 
directly. These are y^ , y^, and in Figures 11, 12, 13, 

and 14, respectively. These figures indicate that the largest 
variations of the metrics occiir near the trailing edge and at 
the far field boundary. They are quite uniform near the leading 
edge. This indicates that the mesh has been sufficiently re- 
fined there to reduce the solid body angle variation to accept- 
able levels. It may be possible, however, that the extra 
resolution is needed to resolve the peak of stagnation density 
adequately. 

Based on these results one would expect most of the mesh 
induced truncation errors to occur near the trailing edge and 
at the far field boundary. That this is indeed the case is con- 
firmed in Figure 15. Shown here is the second group of terms 
of equation (35) in terms of percentage of (the numerically 
computed density) . The mesh induced errors occur mostly near 
the trailing edge and the far field boundary. They are on the 
order of 2 to 3%. 

The solution induced truncation error (first group of terms 
of equation (35)] is shown in Figure 16. These errors are much 
larger. They are on the order of 1-2% near the leading edge, 

5% at the shock wave, <5% at the trailing edge, and 2% at the 
far field. The total error is shown in Figure 17. As can be 
seen some of the errors at the trailing edge and far field have 
been cancelled. This is, however, a fortuitous circumstance and, 
in general, should not be expected to occur. The errors are 
still approximately 2% at the leading and trailing edges and 
5% at the shock wave. 
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One concludes that the mesh of Figure 1 is probably more 
than adequate for controlling the mesh induced errors, if 5% 
solution errors are acceptable. If drag calculations are to be 
attempted with this mesh and numerical procedure, then 2% accu- 
racy at the trailing and leading edge is probably not acceptable. 
The errors would need to be reduced by further mesh refinement. 

RECOMMENDATIONS 

Based on the results of this study the following recommen- 
dations are offered; 

1) Generate a mesh based on acceptable variations of either 

a) x^, y^ , x^, y^ for first order differential systems, or 

b) 9 - 52 ' ^12' second order differential 

systems. 

Acceptable variations can be checked by plots made as shown by 
Figures 4 through 11. 

2) Once the converged numerical solution has been obtained 
{hopefully the scheme is still stable) do the truncation error 
analysis and nimerically compute the mesh and solution induced 
truncation errors. If these are not at acceptable levels, then 
a new mesh needs to be generated until acceptable error levels 
are obtained. (It should be standard practice for any numerical 
computation that the truncation errors are also computed.) 

3) For higher numerical schemes the variations of the mesh 
criteria must be less than for lower order schemes. 

4) For dynamic meshes the variations of the mesh criteria 
must again be less than for stationary schemes. 

EXTENSION TO THREE-DI.MENSIONS 

The extensions of the above results is straight forward in 
the sense that one can readily identify the additional parameters 
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that can make up the mesh criteria. For the 3-D systems, there 
are six components to the metric tensor, i.e., 


ds 


2 


i,j = 1,2,3 

g^^^dC^ + 2g^2‘^Cdn + 2g^3dCdc 
+ 2g23dndC + 922'^’^ * 533 '^' 


and there are nine metric terms, i.e. 


S (x,v,z) 
3 (5,n,C) 



If the metric tensors are used as the mesh criteria, then 
three solid body rotation angles need to be included to completely 
describe the mesh properties. 


I'Thile the additional elements for the 3-D case are obvious, 
deriving the truncation errors is probably impossible to do by 
hand. The use of machine algebra would be mandatory. 


CONCLUSIONS 

In this study an attempt has been made to establish some 
criteria for generating finite difference grids. However, it 
has been shown by truncation error analysis that it is not pos- 
sible to do so independently of the numerical solution obtained 
on the mesh. The study has been limited to mesh criteria which 
reduce the truncation errors. These truncation errors depend on 
the order of the differential system, the solution, the form of 
the transformed equation, the order and type of numerical scheme. 
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and the mesh itself. The results of this study show that prob- 
ably the best strategy is to follow the proced\ire given in the 
recommendations section. The mesh cannot be judged by itself. i 

Its use must also be considered. 
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APPENDIX A 

In this appendix the modified equation for a higher order 
numerical solution to the system of equations (13) is derived. 
The numerican scheme is the leap frog scheme. Although this 
scheme is in general unstable > it is second order accurate both 
spatially and temporarily, and thus is represents a typical 
second order scheme. It is much simpler to analyze than say 
MacCormack's scheme. 

The leap frog scheme is given by 


. „n-l ^ - £? 


2At 


2A5 


zi + h = 


2An 


(A-11 


The modified equation for this scheme and system (13) is 


+ h + ^3^ + ^£35 + ^ 3 n " 


Eliminating the higher order time derivatives obtains, 
after considerable algebra, the final modified equation. 


t . 

Wt + 

® [f + 

3C 3! ^2C 

+ ^ g + 

3n ^ 

TT^2n 

^ ^2w® * ® 

4 

• 


■ 



l» 



+ 


f s • 

^2w n 
n 

s + 2f s- • s + 2f 

n ww^ C ww^ 

® 

1 

1 

i 

+ 

n 





1 

4 

- 

— 35 t[^ 2 w" • 


5 * S5 + 

^2w ®n • ^ "5 


‘ (Continued on next page) 
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•>„‘-=,) ♦ * V‘-^n> Vr‘'=n!’ 


* * \/-^.n> 


] - • 


^ * '>2Wj=C • '{ 


* >'2w % • ^ * ^'’w =5 • = * 2h • s * 2h„ „ 3c • s„ 

n s n V n 


* »w„3 • ns * >’v.v"r."5 • "55 * V 55 ''' ■ * '’"”3." ■ "5n 


C5 


C 55 


5n 


^ \"'5n"" * ^ "■ Vnn'^ * 


+ h s • s 

w^w n nn 

n nn 


= O(A^) 


(A-2) 


where 


and 


s = ^ g^ + h 


■^ = 35 


fs+f s, +f s 
" '"n ^• 


3n 


g s + g Sr + g s_ 

5 n 


♦ V * "v =5 * \ * K.Jn 

t >1 <; 5. 


+ h s-„ + h s„„ 
w. ^n w nn 
5n nn 
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The following observations are made: 

1. The time dependent part makes an awful mess of the 
situation. V7e may want to consider steady solution first. But 
they also show why adaptive meshes can slow convergence to steady 
state and cause other difficulties. 

2. There are several flux Jacobian metrices, i.e.. 
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/ h„ , 
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and the higher order Jacobians. 

The terms involving are due to the solution, those in- 
volving etc. are due the metrics (i.e., mesh). The terms 

involving both (witli the higher order Jacobians) are due to both 
the solution and the metrics. 

3. The terms involving s contain the effects of both the 
solution and the mesh. 

4. The observations made for the first order scheme apply 
here also. 

To sim.plify the modified equation to a reasonable level 
assime that the mesh is fixed with time, but still consider time 
variations of the solution. In this case t.he modified equation 
reduces to 





s 


+ 






s + 


a (- 
-w 


s_) 




(A- 3) 
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where 

s = + h 

and 

-=T = ^[^W • • ®] ^ ‘'w • ® 

The modified equation for the leap frog scheme on a fixed 
mesh given above is similar to that of the first order scheme, 
equation (19) . The truncation errors are of course of one 
higher order. Thus it can be concluded that the higher order 
derivatives of the meshes must be smoother for higher order 
schemes than for lower order schemes. To see this consider the 
steady state form of equation (A-3) . 


^ 3^! ^35 "3T^3n ~ ^ 


(A-4) 


Now need to specialize this modified equation to either the SCL 
or WCL form. Doing the SCL form first. 

SCL 


f - y„f - 
g = -y^f + x^g 
h = 0 . 
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(''n^ - * *55)^ * ^fn3c^ - %35’ * 

^*112555 ^ ’^n?^25 ■ ^*n5^2E * ^n^SC " V3sJ 

* ^-j-ycsn^ "■ *S3n^ - 5yc2n^n * 3«;2n’n ' ^''En‘2n 
- y?^3n * ’‘553„| = "(4^1 



z 


(A-51 


The error terms of this equation should be compared with the 
first order scheme error terms given by equations C22) and C23). 
Similar results are obtained for the WCL form of the transformed 
equations. 

WCL 


w = w 




g 

h 


J~^(-y^f + x^g) 


(J'V) - (J"^y.) 

-(J~^x„) + (J~^x^) 

^ Z ^ n 


The WCL modified equation becomes 




1 / \ 

- x^g)J 

+ 

J" (-yrf ^ 

. 


T* r- 
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where 

J ‘ 

IS ■ 

^rV 

.6 n 


(■“' ^•c)^3n * (■’ ^’'!;)53n| 

(A- 6) 


The error terms of this equation should be compared with 
the first order scheme error terms given by equations (2*1) and 
( 25 ) . 
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Figure 3. Carpet plot of density soVution as solved by TAIRC (ref. 17 ) on lower 
half of NACA 0012 airfoil^ » 0.80/ ct = 0®. Viewpoint from 
near the trailing edge and far field ooundary 
(lower left hand corner of figure 1.) 









Carpet plot of length parameter norma 
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Figure 10. Carpet plot of cell aspect ratio normalize' 





ORIGINAL PASH lt> 
OF POOR QUALITY 



Carpet plot of normalized by 0.1773. 
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